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We calculate analytically the past asymptotic decay rates close to an initial singu- 
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larity in general Go spatially inhomogeneous perfect fluid models with an effective 
supported by numerical simulations in a special class of G2 cosmological models. 



equation of state which is stiff or ultra-stiff (i.e., 7 > 2). These results are then 



Our analysis confirms and extends the BKL conjectures and lends support to recent 
isotropization results in cosmological models of current interest (with 7 > 2). 
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The conjectures of Belinskii, Khalatnikov and Lifshitz (BKL) [1] assert that the structure 
of space-time singularities in general relativity (GR) have the following properties (for a 
perfect fluid with 7 < 2, where (7 — 1) = p/ p is denned as the ratio of the pressure p to 



the energy density p): 1. Each spatial point evolves towards the singularity as if it were a 



spatially homogeneous cosmology. 2. Space-times with non-stiff matter, 7 < 2, have the 
property that asymptotically close to the singularity matter is not dynamically significant 
and the singularities in generic four dimensional space-times are space-like and oscillatory 
(mixmaster behavior). 3. In the case of stiff matter, 7 = 2, the matter is not insignificant 
near the singularity (leading to non-oscillatory behavior) and generically have anisotropic 
singularities which are space-like and non-oscillatory. 

We shall investigate the BKL conjectures in cosmological models with an effective equa- 
tion of state 7 = 2 and extend them to models with 7 > 2 by calculating the past asymp- 
totic decay rates in general (Go) spatially inhomogeneous models. These results are then 
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supported by a numerical analysis of the behavior of spatially inhomogeneous solutions to 
Einstein's equations near an initial singularity in a special class of Abelian G 2 spatially 
inhomogeneous models. 

There are a number of cosmological models of current physical interest which have an 
effective equation of state 7 > 2. Although a complete fundamental theory is not presently 
known the phenomenological consequences can be understood by studying an effective low- 
energy theory, which leads to the introduction of additional fields (e.g., scalar fields) in the 
high curvature regime close the Planck time scale. Scalar fields are believed to be abun- 
dant and pervasive in all fundamental theories of physics applicable in the early Universe, 
particularly in dimensionally reduced higher-dimensional theories . In addition, scalar 
field cosmological models are of great importance in the study of the early Universe, par- 
ticularly in the investigation of inflation 3], 0] and "quintessence" scalar field models 



consistent with observations of the present accelerated cosmic expansion p. Superstring 

theory represents the most promising candidate for a unified theory of the fundamental in- 

rn 

teractions, including gravity |2j. It is widely believed that eleven-dimensional supergravity 
represents the low-energy limit of M-theory Q|. In the low-energy limit, to lowest-order in 
both the string coupling and the inverse string tension, all massive modes in the superstring 
spectrum decouple and only the massless sectors remain, which are determined by the cor- 
responding supergravity actions. A definitive prediction of string theory is the existence of 
a scalar field, known as the dilaton, interpreted as a modulus field parametrizing the radius 
of the eleventh dimension. There are two further massless excitations that are common to 
all five perturbative string theories, namely the metric tensor field (the graviton) and an 
anti-symmetric form field. When the higher-dimensional metric is compactified, additional 
form fields and scalar moduli fields are produced. Thus the low-energy effective action of 
the theory essentially reduces to GR plus massless scalar fields. A massless scalar field (or 
moduli field etc.) close to the initial singularity has an effective equation of state 7 = 2. 

Models in which 7 > 2 arise naturally in ekpyrotic and cyclic cosmological models, 
which have a big crunch/big bang transition with a contraction phase dominated by a scalar 
field with 7 > 2 to the future [10j. In particular, it was shown [11] that if 7 > 2, chaotic 
mixmaster oscillations due to anisotropy and curvature are suppressed and the contraction is 
described by a spatially homogeneous and isotropic evolution. This result was subsequently 
generalized to theories where the scalar field couples to p-forms, and it was also shown that 
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Z 2 orbifold compactification also contributes to suppressing chaotic behavior. Indeed, it was 
concluded that chaos is avoided in contracting heterotic M-theory models if 7 > 2 at the 
crunch. 

There is currently great interest in higher- dimensional gravity theories inspired by string 
theory in which the matter fields are confined to a 3+1-dimensional "brane- world" embedded 
in higher dimensions, while the gravitational field can also propagate in the extra bulk 



dimensions 12]. There has been particular interest in the dynamics of the Universe at early 
times in Randall-Sundrum-type brane- world cosmological models [13]. A unique feature 
of brane cosmology is that p 2 dominates at early times, leading to an effective equation 
of state parameter with 7 > 2, which will give rise to completely different behavior to 
that in GR. The cosmological implications of brane world models have been extensively 
investigated Q]. In particular, it was found that an isotropic singularity is a past-attractor 
in all orthogonal Bianchi models [15]. Moreover, the asymptotic dynamical evolution of 
spatially inhomogeneous brane-world cosmological models close to the initial singularity 
was studied numerically and it was found that there always exists an initial singularity, 
characterized by the fact that spatial derivatives are dynamically negligible, which is isotropic 
for all physical parameter values. The numerical results were supported by a qualitative 
dynamical analysis and a calculation of the past asymptotic decay rates |l6| . 

Andersson and Rendall [17] proved that a generic inhomogeneous cosmology with a stiff 
fluid or a scalar field tends to a velocity-dominated solution, which is the Jacobs solution 
[lsl p. 1426], along individual timelines. For a generic inhomogeneous cosmology with 
7 > 2, we will show that it tends to a flat isotropic solution 1 along individual timelines. 



Note that these solutions are self-similar |2l( . In this paper we shall show that these results 
are supported by the calculations of past asymptotic decay rates in general G models (with 
7 > 2) 2 and numerical simulations in a class of G2 models with one tilt degree of freedom. 



Namely the Binetruy, Deffayet and Langlois solution which is essentially the flat Friedmann-Lemaitre 
solution 0, Ch. 2] with 7 > 2. 
2 We assume that a massless scalar field (or a massive scalar field close to the initial singularity) can be 



modelled by a perfect fluid with a stiff equation of state. This is justified in detail in 



smgr 

:i7]. 
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II. ASYMPTOTIC DYNAMICS AT EARLY TIMES FOR G 2 COSMOLOGIES 

We first consider the dynamics of a class of G 2 spatially inhomogeneous cosmological mod- 
els with one spatial degree of freedom. The governing system of evolution equations consti- 
tute a system of autonomous partial differential equations in two independent variables. We 
follow the formalism of |22j which utilizes area expansion-normalized scale-invariant depen- 
dent variables, and we use the separable area gauge to consider analytically and numerically 
the asymptotic evolution of the class of G 2 cosmologies with one tilt degree of freedom 3 near 
the cosmological initial singularity. For recent works on these models (with 7 < 2), we refer 

early times can be derived by following the analyses in [2J| and [la] by exploiting asymptotic 
silence. 4 The results for G 2 models with one tilt degree of freedom are given below. 

A. Case 7 = 2 

For the case 7 = 2, we assume that the following conditions hold uniformly for open 
intervals of x. 



Cx: lim (E 1 1 ,A,N x ,N^E x ,E 2 ,v) = 0, lim (£+, £_) = (E+(x), £_(*)) 

t^r — OO t— + — OO 

C 2 :d x (E 1 1 ,A,N x ,N.,E x , S2, v, £+, £_) are bounded as t — > — 00. 
C 3 : V = 0(f(t)) implies d x V = 0(f(i)) (asymptotic expansions in time 
can be differentiated with respect to the spatial coordinates). 

The decay rates as t —>■ — 00 are then given by 

{E X \A) = e 2t [{E x \A) + 0(f)] , (£+,£_, fi) = + 0(g) (1) 

E 2 = e M [S 2 + O(g)) , S x = e fc2 *[S x + 0(h)] , AL = e fc3 '[AL + 0(h)) (2) 
N x =e 2t [-E 1 1 d x t_t + N x +0(h)} , v = e 2t [-\E 1 1 d x \nQt + v + 0(h)] (3) 



3 This class is described by the area expansion rate /3 = H — ^cn eq. (2.103)] and the normalized 
variables (Ei 1 , A, N x , N-, S x , S 2 , S+, S_, 17, v) defined with respect to a time-like congruence that corre- 
sponds to the separable area gauge. Ei 1 is the frame coefficient, A, N x , iV_ are the components of spatial 
curvatures, S x , £9, E_i_, E_ are the shears, f2 is the normalized density parameter, and v is the tilt of the 
perfect fluid. See |23j, Appendix D] for the governing equations, where we set A — and evolve j3 instead 
j2.'il £xl (2.108)]. We remind the reader that the logarithmic time variable t = ln£ is used. 

4 See 26] for the notion of the silent boundary and its role in past asymptotic dynamics. 
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where E+ := \[l - E 2 - Cl], and 

fci = -3E+ + \/3E_, fc 2 = -2^3E_, k 3 = 2(1 + \/3E_), (4) 
/ = e 4 * + e 2klt , = te 44 + e 2fel * + e 2M + e 2k3 \ h = t(e 4t + e 2M + e 2k * + e 2fc3 *)- (5) 

The area expansion rate (3 and its associated deceleration parameter satisfy 

P = 6-^-^0 + 0(9)}, q = 2 + 0(f). (6) 

The ten hat variables above are functions of x, and satisfy the two constraints 

= E 1 1 d t hi$ + r , = E 1 1 d x In E 2 - r + 3A-V3N x ] , (7) 

where f = — 3AE + — 3iV x E_ + 3iV_E x — 3vQ. That leaves eight hat variables, the same 
number as when we specify the initial conditions for numerical simulations. 5 To ensure the 
convergence of N_, E x and E 2 , the exponents ki, k 2 and k 3 must be positive for all x. This 
implies that the attractor is confined within the region 

E H >-^(1 + E^), E H <0, E H >v / 3E^ (8) 
v3 

inside the Kasner circle (triangle I in Figure Et) in the state space of Hubble-normalized 
variables, where 

Note that this restriction is gauge-dependent. 

The area expansion-normalized Weyl scalar W is given by 

W 2 = |(E + - E 2 _) 2 + |(1 - 3E + ) 2 E 2 _ + 0{g + t 2 e u ) . (10) 

Figure E shows the hat variables as computed at t = —50 and t = —100 for a numerical 
run (e.g., we plot e -2 '^ 1 , e~ 2 *(iV x +tEi 1 d x 'E-), etc). 6 Since the plots at both times coincide, 



5 Only six of the hat variables are essential, since we can use the remaining temporal gauge freedom to set 
A = 0, and parameterize x to set E\ l = 1. 

6 For numerical simulations, we use CLAWPACK. Since it handles exponential growth rather inaccurately, wc 
evolve the following variables instead: 

ln^ 1 , A/E^, E_, e- 2t N Xl lnS x , lniV_, lnfi, e~ 2t artanh(v), lnS 2 , La/3. 

Compare with j^J Appendix D]. One should ensure that S x and N- are positive during the simulation. 




FIG. 1: The hat variables computed at t = —50 and t = —100 coincide, thus confirming the decay 
rates. 

the decay rates are confirmed. Although the constraints tend to zero, they do not tend to 
zero fast enough and equation (JJJ) is not satisfied to a sufficient degree of accuracy (this 
numerical problem can be solved by standard methods). Nonetheless, the numerical results 
are good enough for confirming the decay rates. The numerical simulations consequently 
confirm the calculations of the decay rates, provide evidence for the conditions C\-Cs, and 
hence lend support to the conjectures formulated above. 

We use 32768 grid points with periodic boundary condition, and run the simulation from t = to t = —100. 
For the initial condition, we set 7 = 2, 0o = 1 and 

e = 0.1, (V) = 1, A) = 0, (£_)„ - -0.2, (E x ) - 0.2, (jV_) = 0.2, Q = 1, (S 2 ) = 0.4 

in conjunction with Q, eq. (9.1)]. 
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B. Case 7 > 2 

For the case 7 > 2, we assume that the following conditions hold uniformly for open 
intervals of x. 



d: lim {E 1 1 ,A,N x ,N_,Z + ,Z_,Z x ,E 2 ,v) = 0, 

C 2 :d x {Ex\A,N x ,N„,Z x , E 2 , t>, S + , E_) are bounded as t — > — 00. 
C3 : V = 0(f(t)) implies d x V = 0(f(t)) (asymptotic expansions in time 
can be differentiated with respect to the spatial coordinates). 

The decay rates as t — > — 00 are then given by 

{Ex 1 , A,N_,N x ,v) = r ] [(Si 1 , A, N„,N x ,v) + O(0] (11) 
(E + , E_, E x , E 2 , - 1) = £ [(£+, E_, E x , E 2 , -2E + ) + O(0] (12) 

where 

^ = e |(37-2)* ; e=e |(7-2)*_ (13) 

The area expansion rate (3 and its associated deceleration parameter satisfy 

{3 = e-ht0 + O (O} , (? = 1(37 -2) +0(0- (14) 

The ten hat variables above are functions of x, and satisfy the two constraints 

= Ex x d x ln/3 - §71) , = S 1 1 <9 ;r lnE 2 -(-f7{) + 3i- v^iVx) . (15) 

That leaves eight hat variables, the same number as when we specify the initial conditions 
for numerical simulations. 

The area expansion-normalized Weyl scalar W is given by 

w 2 = \{±i + e 2 _ + e 2 x + ±De + o(e) ■ (i6) 

Figure 121 shows the hat variables as computed at t — —50 and t = —100 for a numerical 
run (e.g. plot E^/r], etc). 7 That the plots at both times coincide confirms the decay rates. 8 
Therefore, we have presented numerical evidence for isotropization in these models. 



7 For the case 7 > 2, we evolve the variables 

ln^x 1 , A/E^, E_/£, £ x /£, N_/r,, N x / V , (fi - artanh(u)/ry, lnE 2 , ln/3. 

We use the same initial condition as before, except with 7 = 2.1. 

8 We again note that numerically the constraints do not tend to zero sufficiently fast. 
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Ei 1 A £+ x1Q -3 S 2 





< x < 2;t < x < 2ji < x < 2x 



FIG. 2: The hat variables computed at t = —50 and t = —100 coincide, thus confirming the decay 
rates. 

III. ASYMPTOTIC DYNAMICS AT EARLY TIMES FOR G COSMOLOGIES 

It would also be of interest to investigate general inhomogeneous (Go) models with 7 > 2 
close to the singularity. We refer to 26( for the equations in the separable volume gauge, 
using Hubble-normalized variables (where we set A = 0, R a = 0). We give the asymptotic 
decay rates below. We hope to be able to numerically simulate the Go cosmologies in future 
work. 9 



9 Simulation of Go cosmologies is expensive and technically difficult. Simulations of vacuum Go cosmologies 
with 50 grid points for each of x l has been carried out recently [2^ . 
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A. Case 7 = 2 

For the case 7 = 2, we assume that the following conditions hold uniformly for open sets 
of x\ 

Ci : lim (E a \ A a , N a/3 ,v a ) = 0, lim S a/3 = t aP 

C2 '■ d x (E a l , A a , N a p, v a , £ a( g) are bounded as t — > —00. 
C3 : V = 0(f(t)) implies d x V = 0(f(t)) (asymptotic expansions in time 
can be differentiated with respect to the spatial coordinates). 

Without loss of generality, we perform a spatially-dependent rotation to set £ a/ 3 = for 
a ^ f3. The decay rates as t — > —00 are then given by (no summation over repeated indices 
below) 

Ej = e {2 - taa)t [E a l + 0(F)] , r a = e^~ taa)t [f a + 0{tF)\ (17) 
A a = e^-^HEjd^ t + A a + 0(tF)} (18) 
N aa = e^^iN^ + 0(tF)\ (19) 

N af3 = e V-%»)t\_ E j d . e lS{ a j : (3) 8 t + fiafl + QQpy ? ^ a ^ (20) 

fi = fi + C(F), w a = e (2 - f: «« ) *[-|( J B Q ^lnQ-2r a ) t + £ a + O(tF)] (21) 
S Qa = t aa + 0(F) , S a/3 = 0(F) , a ^ (3 (22) 

where 

F = ^2(2-,+ )* + e 2(2+2 S _)t s = max S^, sJf) = min E m . 

+ V ' a=l,2,3 V 7 a=l,2,3 

The Hubble scalar and the deceleration parameter satisfy 

H = e- 3t [H + 0{F)} , g = 2 + 0(F). (23) 

The twenty eight hat variables above 10 are functions of x\ and satisfy the following sixteen 
constraints 

f a = -E a l di In H , = E a l di± aa + (2 - ± aa )f a - 3A a ± aa - eJ^NpgtJ + 6Q£ Q , (24) 
Cl = l-i (Sn 2 + S22 2 + S33 2 ) , = 2{E [ Jd j - r [a - A [a )Erf - e^N^Ej . (25) 

10 Note that t n + t 22 + £33 = 0, and we can write t aa = diag(-2£+, £+ + V3S-, £+ - Vz£_). 
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That leaves twelve hat variables, of which eight are essential, since we can use the remaining 
temporal gauge freedom to set H = const and make a change of coordinates to set three of 
the Ej's. 

Note that S QQ , satisfy |£ aQ ,| < 2 by definition. The exponents of N aa must be positive, 
however, thus requiring that E m > —1; i.e., the Jacobs solutions are restricted within the 
triangle 

S + <^, --^=(l + S + )<S_<i=(l + S + ). (26) 

See triangle I in Figure Eb- 

The Hubble-normalized Weyl scalar W is given by 

W 2 = |(E+ + t\ - E 2 ) 2 + \{l - 2£ + ) 2 £ 2 + 0(tF) . (27) 



The work of Andersson and Rendall |17( is more rigorous, but also assumes more (e.g., 
analyticity) . Our results complement their analysis by providing more details on the spatial 
curvature variables and on the big O terms, and is easier to apply and interpret. We have 
also discussed the restrictions on the Jacobs disk. 

B. Case 7 > 2 

For the case 7 = 2, we assume that the following conditions hold uniformly for open sets 
of x l . 

d : lim (E a \A a ,N 0l p,v a ,Y, a p) = 0, 

t^— 00 

C2 '■ d x (E a \ A a , N a p, v a , S Qj g) are bounded as t — > —00. 
C 3 : V = 0(f(t)) implies d x V = 0(f(t)) (asymptotic expansions in time 
can be differentiated with respect to the spatial coordinates). 

The decay rates as t — ■> — 00 are then given by 

(Ej, r a , A a , N a p, v a ) = n [(Ej, r a , A a , N aP , v a ) + C(0] (28) 

v aP = £$00 + oa 2 )] , n = i-l± aP ±^e + o(e + v 2 ) (29) 

where n and ^ are given in (JT3J). The Hubble scalar and the deceleration parameter satisfy 
H = e~fr[H + 0(£ 2 )] , q = \{Z 1 -2) + 0{e) . (30) 
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The hat variables above are functions of x l , and satisfy the following constraints 
r a = -E a % \nH, = 2r a + 3 7 £ Q , = 2{E [ Jd j - f [a - A [a )E^ - e a/3S N S "> 'Ej . (31) 
The Hubble-normalized Weyl scalar W is given by 

w 2 = § ■ \t aP t^e + o(e) . (32) 

We comment that the results for G 2 and Go cases differ slightly, due to the difference in 
both the temporal and spatial gauges employed. 

IV. DISCUSSION 

We have shown analytically (and confirmed numerically for the G2 models) that for the 
case 7 = 2, a subset of the Jacobs solutions (triangles I in Figure EJ) are locally stable into 
the past, and for the case 7 > 2 the flat Friedmann-Lemaitre solution is locally stable into 
the past. These results confirm and extend the BKL conjectures, and complement the work 
of Andersson and Rendall 17]. They are of also of current physical interest; for example, 
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they lend support to previous isotropization results 

In applications to specific physical theories, it may be necessary to investigate the ro- 
bustness of these results in the presence of certain additional fields. For example, string and 
M-theories were discussed in order to motivate the existence of a scalar field (dilaton) in 
the early universe. However, in these theories there are also p-form fields present, which can 
lead to oscillatory behavior close to the cosmological singularity j^j. 

We comment that the conditions C1-C3 in Section Til Al can sometimes be violated in a 
neighbourhood of isolated timelines x = const, where formation of spiky spatial structures 
are observed numerically. In a preliminary numerical investigation for the stiff G2 case, 
these structures appear to be similar to the spikes in vacuum models [24] and in models 
with 7 < 2 [23]. We conjecture that the limits for E + and £_ along spike timelines reside 
in the following two triangles: (i) for true spikes, the limits reside in the triangle 

-A(i + E £) <E ^<__L(i + E £) ; S ^>v / 3S^, (33) 

(ii) for false spikes, the limits reside in the triangle 

S^<^(l + S^), £ H >0, S H <-v / 3Sr; (34) 
\/3 
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FIG. 3: a) For the stiff G2 case, the Jacobs solutions that act as the past attractor are confined 
within triangle L The limits along the conjectured true and fake spike timelines are confined within 
triangles II and III respectively, b) For the stiff Go case, the Jacobs solutions that act as the past 
attractor are confined within triangle I. The limits along the conjectured true spike timelines are 
confined within triangles II. 

(see triangles II and III in Figure respectively). We conjecture that true spikes are 
physically real, and are reflected by a discontinuous limit for the Weyl scalar. On the other 
hand, false spikes are artifacts of the rotating frame, and are reflected by a continuous limit 
for the Weyl scalar. It is of interest to study these spikes in more detail. 

Similarly, for the stiff Go case we conjecture that spikes occur when the limits of Y> aa are 
discontinuous. When this occurs, the limits are confined to the three triangles 

< 4 + 2S QQ - , ± aa < -1 , (35) 

for all a 7^ (3 (see triangles II in Figure Eb; also compare with Figure 3]). 
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